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ABSTRACT 

Using two-dimensional simulations of non-radiative viscous rotating black hole 
accretion flows, we show that the flows with a ~ 0.1 — 0.3 self-organize to form 
stationary unipolar or bipolar outflows accompanied by global meridional circulations. 
The required energy comes, with efficiency ~ 0.001 — 0.01, from the matter directly 
accreted onto black hole. Observational implications are discussed. 
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1 INTRODUCTION 

The considerable recent interest in models of geometrically 
thick rotating accretion flows has been motivated by abil- 
ity of these models to explain unusual observational proper- 
ties of some accreting black hole candidates. There are two 
major classes of such models, (a) At low accretion rates, 
the optically thin accretion flow is not able to radiate effi- 
ciently and the internal energy stored in the flow is advected 
into the black hole. Models of such flows have been devel- 
oped by Ichimaru (1977), Rees et al. (1982), Narayan & 
Yi (1994), Abramowicz et al. (1995), and others see in re- 
cent reviews by Narayan, Mahadevan & Quataert (1998) 
and Kato, Fukue & Mineshige (1998). These models are 
called advection dominated accretion flows (ADAFs) and 
have been applied for explanation of spectral properties of 
low luminosity high-energy sources, (b) At high accretion 
rates the flow is optically thick. The liberated binding energy 
is converted to radiation which is trapped by the flow and 
advected into the black hole as shown by Katz (1977) and 
Begelman (1978) for spherical accretion. Rotating accretion 
flows of this type are called thick discs, Polish doughnuts, 
or slim discs. They have been developed by Abramowicz, 
Jaroszyhski & Sikora (1978), Jaroszyhski, Abramowicz & 
Paczyhski (1980), Abramowicz et al. (1988) and others. 

Most of the ADAF models have been constructed in a 
one-dimensional approach which restricts properties of the 
solutions by considering only the vertically-averaged purely 
inward motion: the multidimensional character of flow is 
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missing. Narayan & Yi (1995a) pointed out the possible im- 
portance of polar outflows in ADAFs. Analytic study of ac- 
cretion flows with polar outflows was performed by Xu & 
Chen (1997) and Blandford & Begelman (1999) in a self- 
similar approach, and numerical two-dimensional studies 
have recently been carried out by Igumenshchev, Chen & 
Abramowicz (1996), Igumenshchev & Abramowicz (1999, 
hereafter IA99) and Stone, Pringle & Begelman (1999). 
These numerical studies demonstrated that in the case of 
small or moderate viscosity (a 0.1), non-radiative ac- 
cretion flows are convectively unstable and accompanied by 
irregular bipolar outflows. They do not form powerful un- 
bound winds. At large viscosity (0.1 $ a < 1) flows are 
stable and may have (but do not have to have) strong out- 
flows. 

In this Letter we present results from two-dimension 
axisymmetric hydrodynamical simulations of non-radiative 
rotating accretion flows with large viscosity. The flows form 
powerful unipolar or bipolar outflows and could be station- 
ary or nonstationary, depending both on a and on the adia- 
batic index 7, and self-organize into global meridional circu- 
lation cell(s) . The energy required to support the circulation 
is extracted from matter accreted into the black hole with 
efficiency ~ 10~ 3 - 10" 2 . 



2 NUMERICAL METHOD 

We compute models by solving the non-relativistic hydro- 
dynamical equations 



g+pW=0, 

Ad 

p-£ = -vp + pv$ + vn, 



(i) 

(2) 
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Table 1. Parameters of two-dimension models. 



Model 


a 


7 


outflow 


stability 


€ 


A 


0.3 


5/3 


bipolar 


stable 


0.002 


B 


0.1 


5/3 


unipolar 


stable 


0.012 


C 


0.3 


4/3 


unipolar 


quasi-stable 


0.002 


D 


0.1 


4/3 


bipolar 


unstable 


0.006 



de 

P^ = -PVv + Q, (3) 

where p is the density, v is the velocity, $ = —GM/r is the 
Newtonian gravitational potential for a central point mass 
M, e is the specific internal energy, n is the viscous stress 
tensor with all components included, and Q is the dissipation 
function. We adopt the ideal gas equation of state, P = 
(7 — l)pe, and consider only the shear viscosity, assuming 
the o-prescription, 




(4) 



'P/p is the isothermal sound speed, and Q,k = 
yj GM/r 3 is the Keplerian angular velocity. Details of the 
numerical technique used to solve (l)-(3) in axial symmetry 
were discussed by IA99. 

We use a spherical grid N r x Ng = 130 x 50 with 
the inner radius at Ti n = 3r g and the outer radius at 
Tout = 8000r 9 , where r g — 2GM/c 2 is the gravitational 
radius of black hole. The grid extends from to 7t in the po- 
lar direction. To model the relativistic Roche lobe overflow, 
that governs the flow near to the black hole (Abramowicz 
1981), in Newtonian gravity we adopt absorbing boundary 
conditions at r = r\ n together with the condition of the 
derivatives d(vg/r)/dr and d(v$/r) / dr being zero. The lat- 
ter means that there is no viscous energy flux from the inner 
boundary associated with the (rO) and (r(f>) components of 
the shear stress. At r out we apply free outflow boundary 
conditions by interpolating all dynamical variables behind 

7 'out ■ 

In the calculations we assume a source of matter with a 
constant ejection rate. The source is located around = n/2 
in the vicinity of r ou t. Matter is ejected there with angular 
momentum equal to 0.95 times the Keplerian angular mo- 
mentum. Due to the viscous spread, part of the ejected mat- 
ter moves inwards and forms the accretion flow. The other 
part leaves the computation domain freely through the outer 
boundary. We start the computation of a model from an ini- 
tial state, the choice of which is not crucial, and follow the 
evolution until a quasi-stationary flow pattern is established. 
Typically, this takes a few viscous time scales estimated at 

T out ■ 



3 TWO-DIMENSION HYDRO DYNAMICAL 
MODELS 

Four models with various values of a and 7 are listed in 
Table 1, which also lists the type of the outflow, the stability 
and the efficiency e defined by (8). All of the models have 
powerful outflows, launched at radial distances 10 — 100r 9 . 
Models A and B are stationary, and Model C shows a stable 
time-averaged global flow pattern, perturbed from time to 




500 



Figure 1. The flow pattern in Model B. Distributions of density 
p and the vector r 2 pv are shown in the meridional cross section. 
The vertical axis coincides with the axis of rotation. The black 
hole is located in the origin (0,0). The contour lines are spaced 
with A log p = 0.1. Unipolar outflow is clearly seen in the vectors. 
Only a small fraction of the matter circulating in the computation 
domain directly accretes into the black hole. 



time by hot convective bubbles that originate very close to 
r in . Model D is less stationary, with significant convection 
activity which originates in the innermost part of the flow. 

The flow patterns with bipolar outflows in Models A 
and D are quite similar to those discussed by IA99 (their 
Models 1 and 5, respectively), despite differences in a and 
7. The most interesting feature of Models B and C is the 
unipolar outflow. Figure 1 presents the flow pattern for 
Model B. Matter contained within the calculation domain 
forms a global meridional circulation cell with the spatial 
scale ~ Tout- Most of the stream-lines of the flow start at the 
source of matter at r ou t, go inward, approach some minimum 
radius near to the equatorial plane, turn outward, and leave 
the computational domain through the upper hemisphere. 
In Model B we observe one large circulation cell (toroidal in 
three-dimensions), whereas in Model C the large circulation 
cell co-exists with a smaller one of opposite vorticity. The 
circulation is powered by the one-sided outflow generated 
in the vicinity of black hole. The part of the outflow which 
is close to the polar axis is most efficiently accelerated and 
becomes supersonic at a radial distance ~ 1000r 9 . This part 
of the flow contains a small mass fraction and has outward 
directed veloci ties, whi ch are larger than the escape velocity, 
v r > v esc = yj 2GM/r. Obviously, it can escape to infinity, 
even if cooling processes become efficient at large distances 
r i; 1000r 9 . However, most of the outflowing matter forms 
a 'subsonic wind'. The evolution of this wind at large dis- 
tances will be governed by cooling processes, which are not 
considered in our models. 

What drives such powerful outflows and circulations? 
Obviously, the power is extracted from the rotating accretion 
flow with the help of a mechanism which redistributes energy 
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and momentum between different parts of fluid. In our non- 
radiative viscous models only the shear stress can provide 
such a mechanism. Due to the stress, the inner more rapidly 
rotating parts of the accretion flow pass a fraction of kinetic 
(not only rotational) energy to the outer parts. 

The importance of outward energy transport supported 
by convection in geometrically thick accretion discs has been 
independently recognized by Bardeen (1973), Abramowicz 
(1974) and Bisnovatyi-Kogan & Blinnikov (1977), and first 
studied in detail by Paczyriski & Abramowicz (1982) with 
a follow-up by Rozyczka & Muchotrzeb (1982). However, 
these models do not include inward advection. In later works 
(Begelman & Meier 1982; Narayan & Yi 1995a; Honma 1996; 
Kato & Nakamura 1998; Manmoto et al. 2000) advection 
was included, but the convective outward energy flux was 
found to be weak, and was always dominated by the assumed 
purely inward directed advective energy flux. In the 'self- 
similar' models of ADAFs (Narayan & Yi 1994) the viscous 
energy flux, 

2 r 

E visc (r) = —2nr / (v r U rr + v e U r g + v^Ur^) cos Odd, (5) 



is positive, directed outward, and is exactly balanced by the 
inward advection of energy with the corresponding flux 



E adv (r) = 2nr 2 



pv r 



w 



cos 6d8, 



(6) 



v I 

where W is the specific enthalpy. Thus, the total energy flux 



E — E a dv + E v i 



(7) 



is zero everywhere, E(r) = 0. Abramowicz, Lasota & Igu- 
menshchev (1999) demonstrated that E = is an artifact 
of the assumed self-similarity: ADAFs that obey physically 
reasonable boundary conditions have, in general, B ^ 0. 

Our models have a specific geometry of flow which is 
very different from pure equatorial inflow. In Figure 2 we 
show the flow pattern for Model B in the vicinity of the 
inner boundary. In this figure the arrows show the veloc- 
ity directions, and the ellipse is the projection of the radius 
r = Ra at the equatorial plane. Matter, which crosses the 
equatorial plane inside Ra, accretes into the black hole at 
the rate M. Using the analogy of the flow pattern in Fig- 
ure 2 with Bondi-Hoyle accretion, we call Ra the 'accretion 
radius'. In the case of Model B, Ra « 10r 9 . Model C has 
a similar flow pattern, but Ra varies with time around the 
average value « 80r 9 . 

The flow geometry shown in Figure 2 allows a signif- 
icant transport of energy by shear stress from the matter 
accreted into the black hole to the matter outflowing to the 
upper hemisphere. Indeed, the stream lines of matter which 
is eventually inflowing and outflowing are located close to- 
gether until the accreting matter reaches r ~ Ra- During 
this phase, efficient momentum and energy exchange takes 
place. We shall characterize the outward energy transport 
E by the 'accretion efficiency' 



e = E/Mc, 



(8) 



which is different from the standard definition of radiative 
efficiency, 

trad = L/MC 2 , (9) 

where L is the total luminosity of accretion flow. 
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Figure 2. Two-dimensional structure of accretion flow in the 
innermost part of Model B. The region with radius 20r g is shown. 
The axis of symmetry of the picture coincides with the axis of 
rotation. The central shadow circle indicates the location of the 
inner numerical boundary. Arrows point to directions of motion. 
The ellipse is the projection of a circle of 'accretion' radius Ra at 
the equatorial plane. Matter, which crosses the equatorial plane 
inside Ra is absorbed by black hole. Regions with supersonic 
inward velocities are shown in grey. 



Assuming that matter falls freely into the black hole 
inside Ra and that only the binding energy at the orbit r = 
-Ra is liberated, one can obtain an estimate of the maximum 
accretion efficiency, 

1 r„ 



4Ra' 



(10) 



In the case of Models B and C, formula (10) predicts e w 
0.02 and 0.004, respectively. Both values are significantly 
larger than the prediction of e rac j for ADAFs, e ra£ j 10 -4 
(Narayan & Yi 1995b). 

E(r), E a d v (r) and E v i ac (r) for Model B are shown in 
Figure 3. In a steady state, E must be constant. The varia- 
tion in E seen in Figure 3 is about 10%. This can be partially 
explained by a small non-stationarity of the flow, and small 
inaccuracies in our numerical code, which does not exactly 
conserve the total energy. Test calculations have shown that 
this inaccuracy gives an error of less than 5%. From Figure 3, 
one can see that for Model B, e ~ 0.01, which is only a factor 
of 2 smaller than what is predicted by estimate (10). The 
viscous flux E V i 3C (the dashed line in Fig. 3) is always posi- 
tive and dominates the advection flux E a dv (solid line) inside 
the radius ~ 30r 9 . At larger radii the energy is transported 
outward mainly by advection. Note, that E a d v changes sign 
at r ss 5r g . Inside this radius the inward energy advection 
compensates a fraction of the outward viscous energy flux. 

In all models, the time averaged net accretion rate M(r) 
through successive spheres of radius r, is a constant inside 
~ 1000r 9 to good accuracy. However, the rates of mass in- 
flow Mi n and outflow M ou t = Mm — M are both increasing 
functions of r. In Model B Mi n (r) is well approximated by a 
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Figure 3. Radial dependence of the total energy flux E (heavy 
line), the advective energy flux E at j v (solid line) and the viscous 
energy flux E v i 3c (dashed line) in Model B. The energy fluxes are 
defined by eqs. (7), (6) and (5), respectively. 

power law with index /3 « 0.5 in the radial range 10 — 1000r 9 . 
Models A, C and D show a similar power law behaviour for 
Min(r), but with slightly different indices. Such a fast in- 
crease outward of Mi n and M ut indicates the importance 
of global circulation motions in the dynamics of accretion 
flow in the models considered. Only a small part of the mat- 
ter circulated in the computational domain is accreted by 
the black hole. Most of the matter in the case of (quasi) sta- 
ble Models A, B and C, which starts at the source near to 
the outer boundary, has escaped outside the outer boundary 
after one cycle of circulation. This behaviour is in agreement 
with the results of Blandford & Begelman (1999), but unlike 
these authors, we see powerful outflows only in the case of 
large viscosity, a 0.1. We also note that the radial de- 
pendence of Min and M out in our simulations is in good 
agreement with the results of Stone et al. (1999), although 
in thier models accretion flows are dominated by small-scale 
vortices and circulation motions. 



4 DISCUSSION 

Our numerical models can be applied to the two well-studied 
accretion regimes mentioned in the Introduction: optically 
thick and optically thin. The matter in optically thick flows 
is radiation dominated and characterized by 7 = 4/3, so 
Models C and D can be relevant in this case. In optically 
thin flows the 'effective' adiabatic index ranges from about 
3/2 to 5/3, depending on the strength of the magnetic field 
(see Narayan & Yi 1995b) . 

Flow patterns of 'optically thin' models consist of large- 
scale stable subsonic circulation cell(s) (in the r — 9 plane): 
two equatorially symmetric cells in Model A and one cell 
in Model B. Only a small fraction of the outflowing matter 
forms the unbound supersonic unipolar outflow in Model B. 



unipolar 

supersonic 

outflow . , . 

circulation 




Re 

Figure 4. Schematic illustration of the proposed geometry of an 
optically thin rotating accretion flow accompanied by the global 
meridional circulation of matter. 



A fraction e ~ 0.01 of the total energy of the matter accreted 
by the black hole is extracted and remains in the form of ki- 
netic and thermal energy of the circulated matter. Although 
the present numerical models include no energy losses, one 
can speculate about the fate of subsonic outflows when en- 
ergy losses are important, using the recent results obtained 
by Allen, Di Matteo & Fabian (1999) and Di Matteo et 
al. (1999). They have noted that in optically thin accre- 
tion flows coupled with strong outflows, the radiative energy 
losses are dominated by bremsstrahlung. The most efficient 
cooling in such accretion flows occurs at the maximum ra- 
dius where this flow exists. Taking into account this result 
and the results of our numerical simulations we propose the 
following scenario of black hole accretion accompanied by 
global meridional circulation of matter. Subsonic outflows 
originate near to the black hole and are mainly supported by 
the pressure gradient force while the matter radiates weakly. 
At large radial distances, where the dynamical time scale of 
flows becomes comparable with the bremsstrahlung cooling 
time, the matter cools down and its outward motion cannot 
be supported any longer by the pressure gradient. The mat- 
ter turns back and forms the inflowing part of the circulation 
cell(s) . A schematic illustration of such a flow pattern in the 
case of one circulation cell near to the black hole is shown 
in Figure 4. Assuming that all of the thermal energy carried 
by outflows is radiated away, we can roughly estimate the 
radiative efficiency to be e ra d « e ~ 0.01. Note that this sce- 
nario can only be correct if the thermal instability develops 
slowly at the outer region of circulation cell(s), but this is 
what one would expect. 

Basic properties of accretion flows with circulations can 
roughly be estimated. Let R c be the spatial scale of the cir- 
culation cells which contain matter with the average density 
p c . Comparing the bremsstrahlung cooling time for matter 
with the virial temperature and the dynamical time scale 
td = Rc/Vc (where V c = /3Vk is the mean circulation veloc- 
ity at R c , Vk is the Keplerian velocity and the value of the 
parameter /3 ~ a ~ 0.1 is taken from numerical models), we 
can estimate p c oc R~ 2 and consequently, the mass involved 
in circulations, 




and the bremsstrahlung luminosity, 
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Thus, the less massive and more compact object the higher 
luminosity it has. In a steady state, the external mass supply 
to the circulation cells is equal to the mass accretion rate 
M = L/ec 2 . Without the mass supply, the characteristic 
lifetime of the circulation cells is 

Note, that in estimates (11)-(13), we ignore the mass and 
energy losses due to the supersonic unipolar/bipolar out- 
flows, and due to a wind which starts on the 'outer surface' 
of circulation cells and carries out all of the excess angular 
momentum. 

Estimates (11) and (12) agree quite well with the ob- 
served data for the core of the elliptical galaxy M87 (Allen 
et al. 1999). Also, our finding of unipolar outflows from ac- 
creting black holes can be used to explain the one-sided jets 
observed in M87 and other objects. 
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